!-------------------------------------------------------------LICENSE--------------------------------------------------------------!
!                                                                                                                                  !
!The MAP code is written in Fortran language for magnetohydrodynamics (MHD) calculation with the adaptive mesh refinement (AMR)    !
!and Message Passing Interface (MPI) parallelization.                                                                              !
!                                                                                                                                  !
!Copyright (C) 2012                                                                                                                !
!Ronglin Jiang                                                                                                                     !
!rljiang@ssc.net.cn                                                                                                                !
!585 Guoshoujing Road. Pudong, Shanghai, P.R.C. 201203                                                                             !
!                                                                                                                                  !
!This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License         !
!as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version.             !
!                                                                                                                                  !
!This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of    !
!MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.                        !
!                                                                                                                                  !
!You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software     !
!Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.                                                   !
!                                                                                                                                  !
!-------------------------------------------------------------LICENSE--------------------------------------------------------------!

!==================================================================================================================================|
subroutine neigh_reflux (base)
!==================================================================================================================================|

   use parameters
   use block_define
   implicit none

   type(block), pointer, intent(in) :: base

   integer(4) :: i, j, i_cpu, i_cpu_p1, i_neigh, j_neigh, i_level, i_request, i_tmp, j_tmp, sgn
   integer(4) :: n_reflux, i_end, j_end, nx_tmp, ny_tmp, i_n, j_n, i2, j2, i2m1, j2m1
   integer(4) :: i_recv, i_send, cpu_index_x, cpu_index_y
   integer(4) :: n_request, n_neighs_cpu, pos_x, pos_y
   integer(4), dimension(:), allocatable :: request
   integer(4), dimension(:, :), allocatable :: request_status

   real(8) :: dt_dx, dt_dy

   type(block), pointer :: current, neigh

!----------------------------------------------------------------------------------------------------------------------------------|
   n_request = count(n_neighs_reflux_sendrecv_x(:, cpu_id_p1) .gt. 0) +                                                            &
               count(n_neighs_reflux_sendrecv_x(cpu_id_p1, :) .gt. 0) +                                                            &
               count(n_neighs_reflux_sendrecv_y(:, cpu_id_p1) .gt. 0) +                                                            &
               count(n_neighs_reflux_sendrecv_y(cpu_id_p1, :) .gt. 0)
   allocate(request(n_request * 2), request_status(mpi_status_size, n_request * 2))

!----------------------------------------------------------------------------------------------------------------------------------|
!  collect information for sending data
!----------------------------------------------------------------------------------------------------------------------------------|
   cpu_index_x = 1
   cpu_index_y = 1
   do i_cpu = 0, cpu_number - 1

      if (i_cpu .eq. cpu_id) cycle
      i_cpu_p1 = i_cpu + 1

      n_neighs_cpu = n_neighs_reflux_sendrecv_x(i_cpu_p1, cpu_id_p1)
      do i_send = 1, n_neighs_cpu
         i_level = 1
         current => base
         do while (neighs_reflux_info_send_x(i_level * 2 - 1, i_send, i_cpu_p1) .ne. 0)
            pos_x = neighs_reflux_info_send_x(i_level * 2 - 1, i_send, i_cpu_p1)
            pos_y = neighs_reflux_info_send_x(i_level * 2, i_send, i_cpu_p1)
            current => current%child(pos_x, pos_y)%pblock
            i_level = i_level + 1
            if (i_level .gt. n_levels) exit
         enddo

         i_neigh = neighs_reflux_info_send_x(n_d_l + 1, i_send, i_cpu_p1)
         if (i_neigh .eq. 1) i_neigh = 2
         if (i_neigh .eq. -1) i_neigh = 1

         do j = 1, ny_no_bnd
            j_tmp = j + ng
            neighs_reflux_send_x(1, j, cpu_index_x) = current%fxro_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(2, j, cpu_index_x) = current%fxmx_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(3, j, cpu_index_x) = current%fxmy_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(4, j, cpu_index_x) = current%fxmz_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(5, j, cpu_index_x) = current%fxbx_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(6, j, cpu_index_x) = current%fxby_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(7, j, cpu_index_x) = current%fxbz_bnd(j_tmp, i_neigh)
            neighs_reflux_send_x(8, j, cpu_index_x) = current%fxen_bnd(j_tmp, i_neigh)
!            neighs_reflux_send_x(9, j, cpu_index_x) = current%fxpo_bnd(j_tmp, i_neigh)
         enddo
         cpu_index_x = cpu_index_x + 1
      enddo

      n_neighs_cpu = n_neighs_reflux_sendrecv_y(i_cpu_p1, cpu_id_p1)
      do i_send = 1, n_neighs_cpu
         i_level = 1
         current => base
         do while (neighs_reflux_info_send_y(i_level * 2 - 1, i_send, i_cpu_p1) .ne. 0)
            pos_x = neighs_reflux_info_send_y(i_level * 2 - 1, i_send, i_cpu_p1)
            pos_y = neighs_reflux_info_send_y(i_level * 2, i_send, i_cpu_p1)
            current => current%child(pos_x, pos_y)%pblock
            i_level = i_level + 1
            if (i_level .gt. n_levels) exit
         enddo

         j_neigh = neighs_reflux_info_send_y(n_d_l + 1, i_send, i_cpu_p1)
         if (j_neigh .eq. 1) j_neigh = 2
         if (j_neigh .eq. -1) j_neigh = 1

         do i = 1, nx_no_bnd
            i_tmp = i + ng
            neighs_reflux_send_y(1, i, cpu_index_y) = current%fyro_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(2, i, cpu_index_y) = current%fymx_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(3, i, cpu_index_y) = current%fymy_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(4, i, cpu_index_y) = current%fymz_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(5, i, cpu_index_y) = current%fybx_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(6, i, cpu_index_y) = current%fyby_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(7, i, cpu_index_y) = current%fybz_bnd(i_tmp, j_neigh)
            neighs_reflux_send_y(8, i, cpu_index_y) = current%fyen_bnd(i_tmp, j_neigh)
!            neighs_reflux_send_y(9, i, cpu_index_y) = current%fypo_bnd(i_tmp, j_neigh)
         enddo
         cpu_index_y = cpu_index_y + 1
      enddo
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  send data to other cpu
!----------------------------------------------------------------------------------------------------------------------------------|
   i_request = 0
   cpu_index_x = 1
   cpu_index_y = 1
   do i_cpu = 0, cpu_number - 1

      if (cpu_id .eq. i_cpu) cycle
      i_cpu_p1 = i_cpu + 1

      n_neighs_cpu = n_neighs_reflux_sendrecv_x(i_cpu_p1, cpu_id_p1)
      if (n_neighs_cpu .ne. 0) then
         i_request = i_request + 1
         call mpi_isend(neighs_reflux_info_send_x(:, :, i_cpu_p1), (n_d_l + 1) * n_neighs_cpu,                                     &
            mpi_integer4, i_cpu, tag_int + 1, mpi_comm_world, request(i_request), error)
         i_request = i_request + 1
         call mpi_isend(neighs_reflux_send_x(:, :, cpu_index_x), sys * ny_no_bnd * n_neighs_cpu,                                   &
            mpi_real8, i_cpu, tag_var + 1, mpi_comm_world, request(i_request), error)
         cpu_index_x = cpu_index_x + n_neighs_cpu
      endif

      n_neighs_cpu = n_neighs_reflux_sendrecv_y(i_cpu_p1, cpu_id_p1)
      if (n_neighs_cpu .ne. 0) then
         i_request = i_request + 1
         call mpi_isend(neighs_reflux_info_send_y(:, :, i_cpu_p1), (n_d_l + 1) * n_neighs_cpu,                                     &
            mpi_integer4, i_cpu, tag_int + 2, mpi_comm_world, request(i_request), error)
         i_request = i_request + 1
         call mpi_isend(neighs_reflux_send_y(:, :, cpu_index_y), sys * nx_no_bnd * n_neighs_cpu,                                   &
            mpi_real8, i_cpu, tag_var + 2, mpi_comm_world, request(i_request), error)
         cpu_index_y = cpu_index_y + n_neighs_cpu
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  receive data from other cpu
!----------------------------------------------------------------------------------------------------------------------------------|
   cpu_index_x = 1
   cpu_index_y = 1
   do i_cpu = 0, cpu_number - 1
      if (cpu_id .eq. i_cpu) cycle
      i_cpu_p1 = i_cpu + 1

      n_neighs_cpu = n_neighs_reflux_sendrecv_x(cpu_id_p1, i_cpu_p1)
      if (n_neighs_cpu .ne. 0) then
         i_request = i_request + 1
         call mpi_irecv(neighs_reflux_info_recv_x(:, :, i_cpu_p1), (n_d_l + 1) * n_neighs_cpu,                                     &
            mpi_integer4, i_cpu, tag_int + 1, mpi_comm_world, request(i_request), error)
         i_request = i_request + 1
         call mpi_irecv (neighs_reflux_recv_x(:, :, cpu_index_x), sys * ny_no_bnd * n_neighs_cpu,                                  &
            mpi_real8, i_cpu, tag_var + 1, mpi_comm_world, request(i_request), error)
         cpu_index_x = cpu_index_x + n_neighs_cpu
      endif

      n_neighs_cpu = n_neighs_reflux_sendrecv_y(cpu_id_p1, i_cpu_p1)
      if (n_neighs_cpu .ne. 0) then
         i_request = i_request + 1
         call mpi_irecv(neighs_reflux_info_recv_y(:, :, i_cpu_p1), (n_d_l + 1) * n_neighs_cpu,                                     &
            mpi_integer4, i_cpu, tag_int + 2, mpi_comm_world, request(i_request), error)
         i_request = i_request + 1
         call mpi_irecv(neighs_reflux_recv_y(:, :, cpu_index_y), sys * nx_no_bnd * n_neighs_cpu,                                   &
            mpi_real8, i_cpu, tag_var + 2, mpi_comm_world, request(i_request), error)
         cpu_index_y = cpu_index_y + n_neighs_cpu
      endif
   enddo

!----------------------------------------------------------------------------------------------------------------------------------|
!  waiting for mpi communication
!----------------------------------------------------------------------------------------------------------------------------------|
   call mpi_waitall(size(request), request, request_status, error)

!----------------------------------------------------------------------------------------------------------------------------------|
!  update neigh blocks
!----------------------------------------------------------------------------------------------------------------------------------|
   cpu_index_x = 1
   cpu_index_y = 1
   do i_cpu = 0, cpu_number - 1

      if (i_cpu .eq. cpu_id) cycle
      i_cpu_p1 = i_cpu + 1

      n_neighs_cpu = n_neighs_reflux_sendrecv_x(cpu_id_p1, i_cpu_p1)
      do i_recv = 1, n_neighs_cpu
         i_level = 1
         current => base
         do while (neighs_reflux_info_recv_x(i_level * 2 - 1, i_recv, i_cpu_p1) .ne. 0)
            pos_x = neighs_reflux_info_recv_x(i_level * 2 - 1, i_recv, i_cpu_p1)
            pos_y = neighs_reflux_info_recv_x(i_level * 2, i_recv, i_cpu_p1)
            current => current%child(pos_x, pos_y)%pblock
            i_level = i_level + 1
            if (i_level .gt. n_levels) exit
         enddo

!         if (.not. current%var_allocate_flag) then
!            allocate(current%ro(nx_block, ny_block), current%mx(nx_block, ny_block), current%my(nx_block, ny_block),               &
!                     current%mz(nx_block, ny_block), current%bx(nx_block, ny_block), current%by(nx_block, ny_block),               &
!                     current%bz(nx_block, ny_block), current%en(nx_block, ny_block), current%gx(nx_block, ny_block),               &
!                     current%gy(nx_block, ny_block), current%po(nx_block, ny_block), current%az(nx_block, ny_block))
!            allocate(current%fxro_bnd(ny_block, n_neighs), current%fxmx_bnd(ny_block, n_neighs),                                   &
!                     current%fxmy_bnd(ny_block, n_neighs), current%fxmz_bnd(ny_block, n_neighs),                                   &
!                     current%fxbx_bnd(ny_block, n_neighs), current%fxby_bnd(ny_block, n_neighs),                                   &
!                     current%fxbz_bnd(ny_block, n_neighs), current%fxen_bnd(ny_block, n_neighs),                                   &
!                     current%fyro_bnd(nx_block, n_neighs), current%fymx_bnd(nx_block, n_neighs),                                   &
!                     current%fymy_bnd(nx_block, n_neighs), current%fymz_bnd(nx_block, n_neighs),                                   &
!                     current%fybx_bnd(nx_block, n_neighs), current%fyby_bnd(nx_block, n_neighs),                                   &
!                     current%fybz_bnd(nx_block, n_neighs), current%fyen_bnd(nx_block, n_neighs),                                   &
!                     current%fxpo_bnd(ny_block, n_neighs), current%fypo_bnd(nx_block, n_neighs))
!            current%var_allocate_flag = .true.
!            if (ini_value_save_flag .ne. 0) then
!               allocate(current%ro0(nx_block, ny_block), current%mx0(nx_block, ny_block), current%my0(nx_block, ny_block),         &
!                        current%mz0(nx_block, ny_block), current%bx0(nx_block, ny_block), current%by0(nx_block, ny_block),         &
!                        current%bz0(nx_block, ny_block), current%en0(nx_block, ny_block), current%az0(nx_block, ny_block))
!            endif
!         endif

!         i_neigh = neighs_reflux_info_recv_x(n_d_l + 1, i_recv, i_cpu_p1)
!         neigh => block%neigh(i_neigh, 0)%pblock
!         if (i_neigh .eq. 1) i_neigh = 2
!         if (i_neigh .eq. -1) i_neigh = 1

!         do j = 1, ny_no_bnd
!            j_tmp = j + ng
!            current%fxro_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(1, j, cpu_index_x)
!            current%fxmx_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(2, j, cpu_index_x)
!            current%fxmy_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(3, j, cpu_index_x)
!            current%fxmz_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(4, j, cpu_index_x)
!            current%fxbx_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(5, j, cpu_index_x)
!            current%fxby_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(6, j, cpu_index_x)
!            current%fxbz_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(7, j, cpu_index_x)
!            current%fxen_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(8, j, cpu_index_x)
!!            current%fxpo_bnd(j_tmp, i_neigh) = neighs_reflux_recv_x(9, j, cpu_index_x)
!         enddo

!print*, current%p_pos, i_neigh
!pause

         i_neigh = neighs_reflux_info_recv_x(n_d_l + 1, i_recv, i_cpu_p1)
         neigh => current%neigh(i_neigh, 0)%pblock
         if (i_neigh .eq. 1) then
            i_neigh = 1
            n_reflux = ng + 1
            sgn = 1
         endif
         if (i_neigh .eq. -1) then
            i_neigh = 2
            n_reflux = nx_block - ng
            sgn = -1
         endif
         j_end = ny_no_bnd / ratio
         ny_tmp = (current%p_pos(2) - 1) * ny_no_bnd / ratio + ng
         dt_dx = min_dt / neigh%dx

         do j = 1, j_end
            j_n = j + ny_tmp
            j2 = j * ratio
            j2m1 = (j - 1) * ratio + 1
            neigh%ro(n_reflux, j_n) = neigh%ro(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(1, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxro_bnd(j_n, i_neigh))
            neigh%mx(n_reflux, j_n) = neigh%mx(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(2, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxmx_bnd(j_n, i_neigh))
            neigh%my(n_reflux, j_n) = neigh%my(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(3, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxmy_bnd(j_n, i_neigh))
            neigh%mz(n_reflux, j_n) = neigh%mz(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(4, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxmz_bnd(j_n, i_neigh))
            neigh%bx(n_reflux, j_n) = neigh%bx(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(5, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxbx_bnd(j_n, i_neigh))
            neigh%by(n_reflux, j_n) = neigh%by(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(6, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxby_bnd(j_n, i_neigh))
            neigh%bz(n_reflux, j_n) = neigh%bz(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(7, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxbz_bnd(j_n, i_neigh))
            neigh%en(n_reflux, j_n) = neigh%en(n_reflux, j_n) + sgn * dt_dx *                                                      &
               (sum(neighs_reflux_recv_x(8, j2m1 : j2, cpu_index_x)) / ratio - neigh%fxen_bnd(j_n, i_neigh))
!               neigh%po(n_reflux, j_n) = neigh%po(n_reflux, j_n) + dt_dx *                                                     &
!                  (sum(current%fxpo_bnd(j2m1 : j2, i_neigh_n)) / ratio - neigh%fxpo_bnd(j_n, i_neigh))
         enddo

         cpu_index_x = cpu_index_x + 1
      enddo

      n_neighs_cpu = n_neighs_reflux_sendrecv_y(cpu_id_p1, i_cpu_p1)
      do i_recv = 1, n_neighs_cpu
         i_level = 1
         current => base
         do while (neighs_reflux_info_recv_y(i_level * 2 - 1, i_recv, i_cpu_p1) .ne. 0)
            pos_x = neighs_reflux_info_recv_y(i_level * 2 - 1, i_recv, i_cpu_p1)
            pos_y = neighs_reflux_info_recv_y(i_level * 2, i_recv, i_cpu_p1)
            current => current%child(pos_x, pos_y)%pblock
            i_level = i_level + 1
            if (i_level .gt. n_levels) exit
         enddo

!         if (.not. current%var_allocate_flag) then
!            allocate(current%ro(nx_block, ny_block), current%mx(nx_block, ny_block), current%my(nx_block, ny_block),               &
!                     current%mz(nx_block, ny_block), current%bx(nx_block, ny_block), current%by(nx_block, ny_block),               &
!                     current%bz(nx_block, ny_block), current%en(nx_block, ny_block), current%gx(nx_block, ny_block),               &
!                     current%gy(nx_block, ny_block), current%po(nx_block, ny_block), current%az(nx_block, ny_block))
!            allocate(current%fxro_bnd(ny_block, n_neighs), current%fxmx_bnd(ny_block, n_neighs),                                   &
!                     current%fxmy_bnd(ny_block, n_neighs), current%fxmz_bnd(ny_block, n_neighs),                                   &
!                     current%fxbx_bnd(ny_block, n_neighs), current%fxby_bnd(ny_block, n_neighs),                                   &
!                     current%fxbz_bnd(ny_block, n_neighs), current%fxen_bnd(ny_block, n_neighs),                                   &
!                     current%fyro_bnd(nx_block, n_neighs), current%fymx_bnd(nx_block, n_neighs),                                   &
!                     current%fymy_bnd(nx_block, n_neighs), current%fymz_bnd(nx_block, n_neighs),                                   &
!                     current%fybx_bnd(nx_block, n_neighs), current%fyby_bnd(nx_block, n_neighs),                                   &
!                     current%fybz_bnd(nx_block, n_neighs), current%fyen_bnd(nx_block, n_neighs),                                   &
!                     current%fxpo_bnd(ny_block, n_neighs), current%fypo_bnd(nx_block, n_neighs))
!            current%var_allocate_flag = .true.
!            if (ini_value_save_flag .ne. 0) then
!               allocate(current%ro0(nx_block, ny_block), current%mx0(nx_block, ny_block), current%my0(nx_block, ny_block),         &
!                        current%mz0(nx_block, ny_block), current%bx0(nx_block, ny_block), current%by0(nx_block, ny_block),         &
!                        current%bz0(nx_block, ny_block), current%en0(nx_block, ny_block), current%az0(nx_block, ny_block))
!            endif
!         endif

!         j_neigh = neighs_reflux_info_recv_y(n_d_l + 1, i_recv, i_cpu_p1)
!         if (j_neigh .eq. 1) j_neigh = 2
!         if (j_neigh .eq. -1) j_neigh = 1

!         do i = 1, nx_no_bnd
!            i_tmp = i + ng
!            current%fyro_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(1, i, cpu_index_y)
!            current%fymx_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(2, i, cpu_index_y)
!            current%fymy_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(3, i, cpu_index_y)
!            current%fymz_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(4, i, cpu_index_y)
!            current%fybx_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(5, i, cpu_index_y)
!            current%fyby_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(6, i, cpu_index_y)
!            current%fybz_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(7, i, cpu_index_y)
!            current%fyen_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(8, i, cpu_index_y)
!!            current%fypo_bnd(i_tmp, j_neigh) = neighs_reflux_recv_y(9, i, cpu_index_y)
!         enddo

         j_neigh = neighs_reflux_info_recv_y(n_d_l + 1, i_recv, i_cpu_p1)
         neigh => current%neigh(0, j_neigh)%pblock
         if (j_neigh .eq. 1) then
            j_neigh = 1
            n_reflux = ng + 1
            sgn = 1
         endif
         if (j_neigh .eq. -1) then
            j_neigh = 2
            n_reflux = ny_block - ng
            sgn = -1
         endif
         i_end = nx_no_bnd / ratio
         nx_tmp = (current%p_pos(1) - 1) * nx_no_bnd / ratio + ng
         dt_dy = min_dt / neigh%dy
         do i = 1, i_end
            i_n = i + nx_tmp
            i2 = i * ratio
            i2m1 = (i - 1) * ratio + 1
            neigh%ro(i_n, n_reflux) = neigh%ro(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(1, i2m1 : i2, cpu_index_y)) / ratio - neigh%fyro_bnd(i_n, j_neigh))
            neigh%mx(i_n, n_reflux) = neigh%mx(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(2, i2m1 : i2, cpu_index_y)) / ratio - neigh%fymx_bnd(i_n, j_neigh))
            neigh%my(i_n, n_reflux) = neigh%my(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(3, i2m1 : i2, cpu_index_y)) / ratio - neigh%fymy_bnd(i_n, j_neigh))
            neigh%mz(i_n, n_reflux) = neigh%mz(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(4, i2m1 : i2, cpu_index_y)) / ratio - neigh%fymz_bnd(i_n, j_neigh))
            neigh%bx(i_n, n_reflux) = neigh%bx(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(5, i2m1 : i2, cpu_index_y)) / ratio - neigh%fybx_bnd(i_n, j_neigh))
            neigh%by(i_n, n_reflux) = neigh%by(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(6, i2m1 : i2, cpu_index_y)) / ratio - neigh%fyby_bnd(i_n, j_neigh))
            neigh%bz(i_n, n_reflux) = neigh%bz(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(7, i2m1 : i2, cpu_index_y)) / ratio - neigh%fybz_bnd(i_n, j_neigh))
            neigh%en(i_n, n_reflux) = neigh%en(i_n, n_reflux) + sgn * dt_dy *                                                      &
               (sum(neighs_reflux_recv_y(8, i2m1 : i2, cpu_index_y)) / ratio - neigh%fyen_bnd(i_n, j_neigh))
!               neigh%po(i_n, n_reflux) = neigh%po(i_n, n_reflux) + dt_dx *                                                     &
!                  (sum(current%fxpo_bnd(i2m1 : i2, i_neigh_n)) / ratio - neigh%fxpo_bnd(i_n, i_neigh))
         enddo

         cpu_index_y = cpu_index_y + 1
      enddo

   enddo

   deallocate (request, request_status)

!----------------------------------------------------------------------------------------------------------------------------------|
   return
end subroutine neigh_reflux
